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Abstract 

In  this  paper  we  analyze  inviscid  aerodynamic  shape  optimization  problems  gov¬ 
erned  by  the  fuU  potential  and  the  Euler  equations  in  two  and  three  dimensions.  The 
analysis  indicates  that  minimization  of  pressure  dependent  cost  functions  results  in 
Hessians  whose  eigenvalue  distributions  are  identical  for  the  full  potential  and  the  Eu¬ 
ler  equations.  However,  the  optimization  problems  in  two  and  three  dimensions  are 
inherently  different.  While  the  two  dimensional  optimization  problems  are  well-posed, 
the  three  dimensional  ones  are  iU-posed.  Oscillations  in  the  shape  up  to  the  smallest 
scale  allowed  by  the  design  space  can  develop  in  the  direction  perpendicular  to  the  flow, 
implying  that  a  regularization  is  required.  A  natural  choice  of  such  a  regularization  is 
derived.  The  analysis  also  gives  an  estimate  of  the  Hessian’s  condition  number  w^hich 
implies  that  the  problems  at  hand  are  ill-conditioned.  Infinite  dimensional  approxima¬ 
tions  for  the  Hessians  are  constructed  and  preconditioners  for  gradient  based  methods 
are  derived  from  these  approximate  Hessians. 
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1  Introduction 


In  recent  years  there  has  been  a  growing  interest  in  solving  optimization  problems  governed 
by  the  Euler  and  the  Navier  Stokes  equations  [l]-[9].  The  new  interest  in  this  classical  field 
[10,  11]  is  due  to  the  increase  in  computer’s  speed  and  improvements  in  algorithms  for  the 
numerical  solution  of  the  flow  equations. 

The  problem  of  designing  a  three-dimensional  wing  requires  solving  an  optimization  prob¬ 
lem  with  many  design  parameters.  Such  a  problem  may  be  computationally  difficult  depend¬ 
ing  on  the  cost  function’s  level  curves  in  the  vicinity  of  the  minimum.  A  measure  for  the 
level  of  difficulty  is  the  condition  number  of  the  Hessian.  The  eigenvalues  of  the  Hessian 
(which  is  a  symmetric  operator)  are  the  curvatures  of  the  cost  function  in  the  principal  di¬ 
rections.  A  large  deviation  in  the  eigenvalues  means  that  the  cost  function  has  level  curves 
which  are  thin  ellipses.  This  is  well  known  in  the  optimization  literature  to  cause  slowness 
of  convergence  toward  the  minimum  for  gradient-based  methods  [12]. 

Aerodynamic  optimization  problems  are  ill-conditioned  as  noted  in  [6,  9],  and  as  will 
be  shown  in  this  paper.  Therefore  gradient  descent  methods  will  be  extremely  inefficient 
especially  when  the  number  of  design  variables  is  large.  A  standard  method  to  overcome  this 
difficulty  is  the  Newton  method  where  the  Hessian  is  computed  explicitly  [12].  The  Newton 
search  direction  is  the  gradient  multiplied  by  the  Hessian’s  inverse,  a  computation  which  is 
impractical  in  aerodynamic  optimization  problems  since  it  involves  numerous  solutions  of 
the  flow  PDEs.  On  the  other  hand  using  low  rank  quasi-Newton  methods,  such  as  BEGS 
[3,  4],  will  result  in  a  deteriorate  convergence  as  the  number  of  design  variables  increases. 
Thus,  a  new  method  is  required. 

Another  difficulty  in  inviscid  aerodynamic  optimization  problems  is  the  ill-posedness 
of  three  dimensional  problems  which  shows  up  as  small  scale  oscillations  in  the  shape  in 
the  direction  perpendicular  to  the  flow.  Such  oscillations  were  observed  in  applications  as 
reported  in  [9].  One  way  to  avoid  these  oscillations  is  to  apply  smooth  finite  dimensional 
representation  of  the  shape  in  the  spanwise  direction.  Another  approach  is  regularization 
by  introduction  of  a  penalty  to  the  cost  function  for  oscillations  in  that  direction.  The 
need  for  penalizing  the  cost  function  in  order  to  remove  oscillations  was  observed  by  [3] 
However,  in  that  case  the  oscillations  were  a  result  of  the  discretization  and  had  no  differential 
counterpart.  Penalization  was  used  also  in  two  dimensions  where  the  differential  optimization 
problem  is  well-posed. 

In  this  paper  we  develop  a  new  approach  to  approximate  the  Hessian  and  its  inverse  for 
optimization  problems  governed  by  PDEs.  Hessian  symbols  were  previously  computed  for 
smoothing  predictions  in  the  development  of  multigrid  one-shot  methods  [13]-[16].  Here,  a 
similar  analysis  is  applied  to  inviscid  flow  problems  including  the  full  potential  and  Euler 
equations  in  two  and  in  three  space  dimensions.  In  Sec.2  the  optimization  problem  is  defined 
together  with  its  small  disturbance  approximation.  The  necessary  conditions  for  a  minimum 
and  their  relation  with  the  Hessian  are  discussed  also  in  the  finite  dimensional  design  space. 
In  Sec. 3  local  mode  analysis  is  presented  to  approximate  the  Hessian’s  symbol.  The  analysis 
is  local  and  involves  freezing  the  coefficients  to  obtain  a  problem  in  half  space  with  constant 
coefficients,  where  Fourier  techniques  are  employed.  In  Sec.  4  the  analysis  is  applied  to  an 
optimal  shape  design  problem  governed  by  the  full  potential  equation.  In  Sec.  5  the  analysis 
is  applied  to  the  Euler  equations,  and  a  symbol  is  obtained  identical  to  the  full  potential 
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case.  In  Sec.  6  the  Hessian’s  symbol  is  analyzed  and  conclusions  are  made  concerning  the  ill- 
conditioning  and  ill-posedness  of  the  problems  at  hand.  In  Sec.  7  regularization  is  discussed 
to  avoid  the  ill-posedness  which  exists  in  three  dimensions  in  the  spanwise  direction.  In  Sec. 8 
preconditioners  are  developed  for  the  small  disturbance  and  the  optimal  shape  problems  in 
subsonic  and  supersonic  flow.  Finally  in  Sec.  9  discussion  and  concluding  remarks  are  made. 


2  The  Optimization  Problem 

A  typical  inviscid  aerodynamic  optimal  shape  design  problem  aims  at  finding  the  shape  of 
a  surface,  e.g.,  airfoil  or  wing,  such  that  the  resulting  pressure  distribution  on  that  surface 
will  minimize  the  least  squares  distance  from  a  prescribed  pressure  distribution.  Let  be 
a  domain  in  and  r(.T)  a  parametric  representation  of  the  part  of  the  boundary  dO.  to 
be  designed.  The  optimal  shape  problem  is  to  compute  the  boundary  position,  T,  that 
minimizes  a  cost  function  defined  on  F,  e.g., 

(2.1) 

where  /  is  a  prescribed  function  and  U  is  the  solution  of  a  PDE  defined  on  fl, 

C{T,U)  =  0. 

2.1  The  Small  Disturbance  Approximation 

For  the  analysis  of  the  Hessian  it  is  enough  to  consider  small  perturbations  of  the  boundary 
r.  In  order  to  further  simplify  the  derivation  we  consider  a  localization  of  the  problem  in  a 
vicinity  of  a  boundary  point  and  study  the  resulting  half  space  problem.  Let  us  introduce 
the  following  notation 

iR”  =  :  r  G  <  a:„}  .  . 

dJR!\.  =  {(,r,r„)  :  x  G  =  0}-  ^ 

We  consider  perturbations  of  the  form 

r  =  F*  -|-  ectn  (2.3) 

U  U*  +  eU  +  0{e'^)  (2.4) 

where  F*  and  U*  are  the  optimal  boundary  shape  and  state  solution  respectively;  n  is  the 
outward  normal  and  £  is  a  sufficiently  small  positive  number.  The  resulting  optimal  control 
problem  is  obtained  using  a  Taylor  expansion  in  e  and  is  discussed  in  more  details  in  the 
next  sections.  The  new  control  problem,  for  the  control  variable  a,  is  defined  in  half  space 
and  is  referred  as  the  “small  disturbance  minimization  problem” ,  namely, 

rmn  F{ct,tj)  (2.5) 

a 

subject  to  the  following  state  equation: 

hil  =  0  on 
BU  =  Cd  on  d]R\ 
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(2.6) 

(2.7) 


and  additional  homogeneous  conditions  on  the  solution  at  infinity.  For  this  control  problem 
the  minimum  is  attained  for  d  =  0. 

The  necessary  conditions  for  a  minimum  of  the  small  disturbance  problem  are  given 
by  the  state  equation  Eq. (2.6-2. 7)  and  two  additional  equations:  design  and  costate.  The 
design  equation  is  a  PDE  defined  on  the  boundary  and  is  denoted  by, 

.A(tf,A,d)  =  0  ondR\,  (2.8) 

where  U  is  the  solution  of  the  state  equation  (2.6-2.7)  and  A  is  the  solution  of  the  costate- 
equation,  i.e., 

LA  =  0  an  R% 

BA  =  Ctf  on  dR\ 

and  additional  homogeneous  conditions  at  infinity. 

We  assume  that  the  perturbation  d  is  composed  of  high  frequencies  and  compute  the 
resulting  solutions,  U  and  A,  in  the  vicinity  of  some  arbitrary  point  on  the  boundary 
The  solution  there  is  approximated  by  a  constant  coefficient  problem  defined  in  half  space 
where  Fourier  analysis  can  be  applied. 

It  can  be  shown  that  for  feasible  solutions  of  the  state  and  costate  equations,  U  =  U{a) 
and  A  =  A(  17(d)),  the  design  equation  residuals  are  equal  to  the  gradient  of  the  cost  function 
with  respect  to  the  design  variables: 

5(d,f7(d))  =  ^(^(d),A(t/(d)),d). 

In  the  vicinity  of  the  minimum  the  following  relation  holds 

g{a,  t7(d))  =  5^(0, 0)  +  Ha  -|-  h.o.t. 

where  we  denote  hy  H  the  Hessian,  i.e., 

H  =  V55(0,0). 

A  Taylor  expansion  of  the  right  hand  side  in  (2.11)  and  comparison  with  (2.12)  results  in 

H  =  AjjUa  +  A^^iJUa  +  Aa.  (2.13) 

The  dimension  of  the  Hessian,  H,  is  determined  by  the  dimension  of  the  design  variable,  d. 
If  d  belong  to  a  finite  dimensional  space  of  dimension  N  then  the  gradient  is  a  vector  of  size 
N  and  the  Hessian  is  a  N  x  N  matrix.  If  d  belongs  to  an  infinite  dimensional  space  (e.g., 
some  function  space)  then  the  gradient  is  an  element  in  an  infinite  dimensional  space  and 
the  Hessian  is  an  infinite  dimensional  operator.  The  Newton  step,  d,  satisfies 

Ha  =  -g  (2.14) 

where  g  is  the  gradient  at  the  given  iteration. 


(2.11) 

(2.12) 


(2.9) 

(2.10) 
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2.2  Finite  Dimensional  Design  Space 

In  applications  it  is  a  common  practice  to  restrict  the  design  space  to  a  finite  dimensional 
subspace  where  the  shape  F,  or  in  the  small  disturbance  model  a,  is  given  as 

a{s)  =  ^a,fj{s),  (2.15) 

i=i 

where  N  is  the  dimension  of  the  design  space,  aj  are  the  design  variables  and  fj  are  fixed 
basis  functions.  The  gradient  for  the  finite  dimensional  case,  g  =  ■  ■  ,gN)-,  is  obtained 

by  a  projection  of  the  infinite  dimensional  one  onto  on  the  finite  dimensional  space,  i.e., 

Jg  =  b  (2.16) 

where  J  is  a  matrix  defined  by 

=  =  (2.17) 

and  the  right  hand  side,  6,  is  a  vector  whose  elements  are  given  by 

h{s)  =  {g{s),fk{s))L2  k  = 

Let  hjk  be  an  element  in  the  matrix  representing  the  Hessian  using  the  basis  functions  ff. 

hjk  =  {HfjJk)L2.  (2.18) 

A  Newton  step  for  minimizing  the  cost  function  consists  of  moving  in  the  direction  d  = 
(di,---,dA7)  given  by: 

h&  =  -b.  (2.19) 

Note  that  this  equation  is  obtained  by  projecting  the  general  equation  (2.14)  onto  a  finite 
dimensional  space  spanned  by 

3  Local  Mode  Analysis 

A  local  analysis  of  the  operator  H  is  done  next.  The  Fourier  analysis  consists  of  analyzing 
the  solution  of  the  following  system  of  equations  in  half  space: 


to 

II 

on 

(3.1) 

BU  =  Ca 

on  d]R\ 

(3.2) 

to 

II 

ii-:i 

on  JR!\. 

(3.3) 

BA  =  CU 

on  dlR'l 

(3.4) 

A{U,A,a)  = 

0  on  dMl, 

(3.5) 
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where  we  replaced  U  by  U  to  account  for  cases  in  which  L  stands  for  a  system  of  PDEs. 

The  computation  of  the  symbol  of  the  Hessian,  (2.13),  near  the  minimum  is  done  by 
considering  a  perturbation  in  the  design  variable  of  the  form 

q(.t)  =  a{k)e‘^'^  (3-6) 

and  calculating  the  corresponding  term  Ha.  The  small  disturbance  solution  U  can  be  rep¬ 
resented  as 

d{x,  x„)  =  't  (3-7) 

J=1 

where  q  equals  the  number  of  boundary  conditions  in  (3.2).  Each  of  the  terms  in  the 
expression  for  U  satisfies  the  equation 

=  0,  (3.8) 

which  implies 

L{lki)^,{k)^l  (3.9) 

where  L{k,  kl)  is  the  symbol  of  L.  Moreover,  Eq.(3.9)  implies  that  for  j  =  1  •  •  •  9, 

detT(^,fc^)  =  0  (3.10) 

and  that  Vj{k)  is  an  eigenvector  of  L{k,  kl)  with  a  zero  eigenvalue.  Substituting  the  expres¬ 
sion  for  U  into  the  boundary  condition  (3.2)  results,  in  the  following  equation, 

j2Mk)Bik,ki)V,{k)  =  C{k)a{k),  (3.11) 

i=l 

where  B{k,  kl)  is  the  symbol  of  B,  and  C{k)  is  the  symbol  of  C.  Introducing  the  matrix 

Wik)  =  {Bil  kl)P,{k), . . . ,  B{1  kl)P,{k)) 


and  the  vector 

Eq.(3.11)  can  be  written  as 

W{k)'fi{k)  =  C{k)&{k), 

(3.12) 

or,  equivalently,  as 

II 

1 

o> 

(3.13) 

where  Uik)  is  defined  by 

(J{x,0)  =  P{k)kk)e^'"  =  ^(fc)e''-" 
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and  where  V'(^)  is  the  matrix  V  —  (V’i(^), . . . ,  V g{k)). 


The  adjoint  equations  are  treated  in  a  similar  way.  The  solution  A  is  represented  by 

A(£  xn)  =  i:  (3.14) 

j-1 

where  q  +  q  =  N  and  N  is  the  degree  of  the  polynomial  in  (3.10).  Let  be  the  roots  of 

det  L{k,kl)  =  0  (3-l'5) 

and  Ej  the  eigenvectors  of  L(^,  k^}  corresponding  to  a  zero  eigenvalue.  The  costate  boundary 
condition  (3.4)  implies 


■tj}i{k)B{k,ki)k,{k)  =  C(k)0(k). 

i=i 

(3.16) 

Introducing  the  matrix 

W(k)  = 

and  the  vector 

m  =  {Mk),...jg{k)) 

Eq.(3.16)  can  be  written  as 

W{k)l3{k)  =  C{k)U{k), 

(3.17) 

or,  equivalently,  as 

A{k)  =  E{k)\y\k)d{k)U{k), 

(3.18) 

where  E{k)  is  the  matrix  :z.{k)  =  (.=.i(fe), . . .  ,.=,^(fc)). 


Substitution  of  (3.13)  and  (3.18)  in  the  symbol  of  the  Hessian’s  expression  (2.13)  results 
in  the  following  formula  for  the  Hessian’s  symbol: 

Hik)  =  AcVW-'^C  +  AKtW~^  CVW-^C  +  Aa.  (3.19) 

4  The  Full  Potential  Equation 

In  this  section  we  apply  the  ideas  discussed  above  to  minimization  problems  governed  by  the 
full  potential  equation.  We  consider  the  problem  of  optimizing  the  shape  of  an  aerodynamics 
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configuration  so  that  the  shape  minimizes  the  deviation  of  the  pressure  from  a  target  pressure 
distribution,  namely, 

mmF{T,p)  =  J^{p  -  p*yda  (4.1) 


subject  to 


>C(T,  <j))  =  div{pV(j))  =  0  on  0 
(l>n  =0  on  r 


(4.2) 

(4.3) 


and  an  additional  boundary  condition  at  the  far-field  (the  notation  (f)n  stands  for  the  outer 
normal  derivative  of  <{)  on  the  boundary). 

The  density,  p,  is  given  by  the  isentropic  density  law  [17]: 


_P_ 

Poo 


1  - 


{V4,y 


2H^ 


1 


(4.4) 


with  7  the  specific  heat  ratio  and  poo  and  Hoc  the  values  of  the  density  and  total  enthalpy 
at  infinity. 

The  pressure,  p,  is  related  to  the  density  p  and  the  speed  of  sound,  c,  by 


P  = 


(4.5) 


where  for  perfect  gases  the  speed  of  sound  is  related  to  the  potential  through  the  relation 


n2  _ 


(7  -  1) 


(Wl 

2  J' 


(4.6) 


The  Mach  number  M  is  given  by 


|V<^|2 


c 


2 


4.1  The  Small  Disturbance  Minimization  Problem 

The  derivation  of  the  small  disturbance  minimization  problem  follows  the  argument  and 
notation  of  Sec.2.1.  We  perform  localization  and  set  the  local  coordinate  system  on  the 
boundary  such  that  the  flow  is  in  the  x-direction  {V4>  -  On  the  perturbed 

boundary  a  Taylor  expansion  gives  [1] 

=p  +  p  +  ap„  +  /i.o.t. 

and 

da’^erturbed  ^  ^  ^ 

R 

where  ^>216  i?i,  R2  are  the  principal  curvatures. 

The  small  disturbance  minimization  problem  is  given  by 

mmF{6:,U)=  f  {{p  -  p*)  +  p{^)  +  apnfdxdy  -  [  ^{p  -  p*fdxdy.  (4.7) 

a  JdR^^  It 
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subject  to 


(4.8) 


(1  —  M'^)^xx  4"  (f>yy  +  ^zz  —  0  IR^. 

The  small  disturbance  boundary  condition  is  given  by 

—  —  (Px&x  4  =  0  (4-9) 

where  the  unperturbed  flow  is  in  the  x-direction.  Since  we  are  analyzing  the  Hessian  for  the 
high  frequencies  it  is  enough  to  consider 

-  ^z  =  4>x^x-,  (4-10) 

The  relation  between  p  and  ^  in  the  cost  function  (4.7)  is  obtained  from  the  relations  (4.4- 
4.6).  We  obtain  that  the  relevant  cost  function  in  terms  of  ^  to  be  minimized  is 

F{a,^)  =  {{p-p*)  -  p<t>xkydxdy  -  a[^{p  -  p*y  -  2pnip  -  p*)]dxdy .  (4.11) 

4.2  The  Adjoint  Equations 

By  standard  variational  calculus  it  can  be  shown  that  the  gradient  of  the  cost  function  is 
given  by 

9  =  <t>xK  -  4  2p„(p  -  p*)  (4.12) 

where  A  is  the  solution  of 

(1  —  M^)\xx  “1"  ^yy  ^ZZ  ~  0  (4.13) 

with  the  Avail  boundary  condition 

-  Aj:  -  2[p<3ix((p  -  P*)  -  pd’x^r)]^  =  0  on  dlR%  (4-14) 

and  6  is  the  solution  of  (4.8)  with  the  boundary  condition  (4.10).  We  also  require  that  the 
solution  is  bounded  for  the  subsonic  case  (M  <  1)  and  that  in  the  supersonic  case  (M  >  1) 
no  waves  propagate  in  the  direction  of  —V(f>{xo),  i.e.,  not  propagating  in  the  negative  x- 
direction.  This  requirement  is  done  so  that  the  solution  of  the  small  disturbance  problem 
will  be  consistent  with  the  far-field  boundary  conditions  of  the  unperturbed  problem. 

4.3  Local  Mode  Analysis 

We  now  go  through  the  analysis  in  Sec..3  in  order  to  compute  the  symbol  H{k)  (see  Eq.3.19) 
using  the  full  potential  state  equations.  Following  a  perturbation  a 

a{x)  =  a{k)e^^'^ 

the  small  disturbance  solution,  can  be  represented  as 

4>{x)  = 

8 


(4.15) 


and,  similarly,  the  adjoint  variable. 


A(f)  =  (4.16) 

Substitution  of  (4.15)  in  (4.8)  and  of  (4.16)  in  (4.13)  results  in  the  following  relations  for 
^3  and  for  k^: 

(1  -  M^)kl  +  kl  +  kl  =  ^  (4.17) 

{I  -  M^)kl  +  kl  +  kl  =  Q.  (4.18) 


These  are  the  analogs  of  equations  (3.10)  and  (3.15)  respectively. 
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Figure  1:  The  supersonic  flow  can  be  decomposed  into  two  waves:  incoming  and  outgoing 
the  plane 

The  choice  of  sign  for  k^  and  kz  is  done  as  follows.  Since  the  half  space  problem  is  related 
to  localization  of  the  original  problem  around  some  point  xq  €  T,  the  solutions  that  we 
construct  in  half  space  should  be  compatible  with  the  far  field  boundary  conditions  of  the 
original  problem.  In  the  supersonic  regime  the  solution  can  be  decomposed  into  two  waves: 
incoming  and  outgoing  of  the  plane  dR%  (see  Fig.l).  In  terms  of  the  local  coordinate  system, 
the  incoming  characteristic  has  a  component  in  the  negative  y-direction  while  the  outgoing 
has  a  component  in  the  positive  y-direction  (both  propagate  in  the  positive  x-direction). 
Since  perturbations  of  the  shape  F  can  not  change  the  far  field  inflow  data,  the  change  in  the 
incoming  characteristic  should  vanish.  Thus,  the  part  of  ^  which  propagates  to  the  negative 
x-direction  is  set  to  zero.  This  implies  that  kz  is  of  opposite  sign  to  that  of  ki  when  kz  is 
real  valued  (i.e.  supersonic  flow).  The  adjoint  variable  has  characteristics  in  the  opposite 
direction  and  therefore  we  require  that  the  part  of  the  solution  for  A  which  propagates  to 
the  positive  x-direction  should  be  set  to  zero.  As  a  result  the  sign  of  kz  is  the  opposite  of 
that  of  kz  if  they  are  real  valued.  In  the  subsonic  regime  kz  and  kz  are  imaginary  therefore 
the  proper  sign  is  positive  for  both  so  that  exponentially  decaying  solutions  are  obtained. 
Therefore  kzkz  =  —  for  both  subsonic  and  supersonic  flow. 

From  the  boundary  condition  (4.10)  we  obtain  a  relation  between  a{k)  and  4>{k),  (as  in 
(3.13)) 

^{k)  =  (4.19) 

kz 
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and  from  the  boundary  condition  (4.14)  (as  in  (3.18)) 

\(k)  =  2,p^4'lrkh  (‘1-20) 

Hence,  the  gradient  (4.12)  in  Fourier  space  is 

g{k)  =  <})xiki\(k)  +  l.o.t.  (4-21) 

Substitution  of  (4.19)  in  (4.20)  and  of  (4.20)  in  (4-21)  results  in  the  relation 

g{t)  =  (4.22) 

and  thus  by  relation  (4.17)  we  obtain  the  following  formula  for  the  symbol  of  the  Hessian: 

H{k)  =  +  (4-23) 

5  The  Euler  Equations 

As  in  the  full  potential  case  we  consider  the  problem  of  optimizing  the  shape  of  an  aerody¬ 
namic  configuration,  subject  to  the  Euler  equations,  such  that  the  optimal  shape  minimizes 
the  deviation  of  the  pressure  from  a  target  pressure  distribution  (4.1).  We  perform  the  anal¬ 
ysis  away  from  shocks  so  that  it  can  be  done  using  a  non-conservative  formulation.  The 
Euler  equations  in  quasi-linear  non-conservative  form  are  given  by 

/  Q  pdx  pdy  pdz  0 

0  Q  0  0  4^ 

£(r,17)=  0  0  Q  0  Idy 

0  0  0  Q  jd, 

\  0  pc'^dx  pc^dy  pc^dz  Q 

where  Q  =  u  ■  V  {u  =  {u.v,w)  denotes  the  velocity  vector),  with  the  solid  wall  boundary 
condition 

u  ■  n  =  0.  (5-2) 

Additional  conditions  that  are  given  at  the  inflow  and  outflow  boundaries  in  terms  of  char¬ 
acteristic  variables  are  not  used  explicitly  in  the  derivation  of  the  approximate  Hessian. 

5.1  The  Small  Disturbance  Minimization  Problem 

Following  the  same  argument  given  in  Sec.4,  the  small  disturbance  cost  function  is  given 
in  Eq.(4.7).  The  state  equations  (5.1)  are  perturbed  in  the  vicinity  of  the  minimum  and 
the  perturbation  variables  solve  the  linearized  Euler  equations  which  up  to  low  order  terms 
are  given  by  the  same  matrix  operator  as  above.  Following  the  localization  and  half  space 
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approximation  we  set  the  local  coordinate  system  on  the  boundary  such  that  the  flow  is  in 
the  x-direction  {u  =  (w,0,0)).  The  small  disturbance  wall  boundary  condition  is  given  by 
(the  perturbation  variables  are  denoted  by  p,  u,  and  p) 

—  w  =  ua^  (b-3) 

where  we  have  omitted  as  before  the  zero  order  terms  in  a. 

5.2  The  Adjoint  Equations 

By  standard  variational  calculus  the  gradient  of  the  cost  function  is  given  by 


g  =  -up{\r  +  c^As)^  -  +  2p„(p  -  p*)  (5.4) 

where  (Ai,  A,  A5)  is  the  solution  of  the  following  system  of  equations  in 

div{Xiu)  =  0  (5.5) 

grad{pXi)  +  grad{u  ■  A)  +  grad{pXs)  =  0  (5-6) 

div{—X)  +  div{X^u)  =  0.  (5-7) 

P 

The  boundary  condition  on  dlR\  is  given  by 

--A4  +  2p(p-p*)  =  0.  (5.8) 

P 


An  additional  requirement  is  consistency  with  the  far-field  boundary  conditions  of  the  un¬ 
perturbed  problem;  i.e.,  the  solution  is  bounded  for  the  subsonic  case  [M  <  1)  and  that  no 
waves  propagate  in  the  direction  of  —u{xq)  in  the  supersonic  case  (Af  >  1). 

Neglecting  zero  order  terms  in  Eqs.(5.5-5.7)  we  arrive  at  the  following  quasi-linear  form: 


Ml  A 

/  g  0  0  0  0  A 

Ai  A 

A2 

0 

0 

A2 

L 

A3 

=  — 

pdy  Q  Q  ^  pc^dy 

A3 

A4 

pdz  0  0  Q  pc^dz 

A4 

As  } 

\  ®  ^  } 

\  ^5  / 

5.3  Local  Mode  Analysis 

Following  the  same  procedure  as  in  the  full  potential  case  we  consider  a  perturbation  in  the 
design  variable  of  the  form 

d(x)  =  d(^)e*'^-"  (5.9) 

with  fc  =  (fci,  ^2)  and  x  =  {x^y).  As  a  result  the  states  and  costates  are  perturbed  by 

U{x)  = 

A(f)  = 
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(5.10) 

(5.11) 


with  U  =  ip,u,p)  and  A  =  (Aj,  A,  A5).  Introducing  the  notation 

Q{k)  =  u  •  k  ==  icki 


and 


L{k,  k-s)  =  i 


(Q 

0 

0 
0 


pki 

Q 

0 

0 


pk2 

0 

Q 

0 


pk3 

0 

0 

Q 


0  \ 


b. 

p 

hz 

p 

hz 

p 


\  0  pc^ki  pc^k2  pc^ks  Q  ) 

the  form  of  the  above  solutions  for  the  perturbation  variables  imply  that 

det  Lik^k^)  =  0 
det  L{k,  ks)  =  0 


(5.12) 

(5.13) 


where  L  is  the  adjoint  of  L.  The  two  relations  result  in  a  fifth-order  polynomial  equation  for 
the  wave  numbers  {ki,k2,k3), 

kr\^(kr\^  -  c\k\  +  kl  -h  kl))  =  0 
ki^u^{ki'^v?  -  c^{kl  -t-  +  ^3))  =  0, 

with  the  roots 

kl  —  k^  = —{1  —  M^)kl  —  kl  and  fci  =  0,  (5.14) 

where  the  Mach  number  is  given  by  Note  that  the  roots  for  ^3  and  for  k^  are 

identical  to  those  obtained  in  the  full  potential  case. 

Let 


hik)  = 

(kipuc  ^^—ki,—k2T—k3,kipu^ 

V2ik)  = 

{kipuc~'^,  —^1,  —^2,  ^3,  k\pv}j 

h{k)  = 

(1,0, 0,0,0) 

i^4ik)  = 

(0, 1,0,0, 0) 

^-.{k)  = 

(0,0, -A;3,  ^12,0) 

kik)  = 

^0,  —c^pki,—c^k2p,  —c^pkz,  kiu^ 

Uk)  = 

(^0,  —c^pki,-c^k2p^  c^ph,  kiuj 

kik)  = 

(-c2,0,0,0,l) 

Uk)  = 

(0,1, 0,0,0) 

Uk)  = 

(0, 0,  —ks^  A:2,  0). 

(5.15) 


(5.16) 


12 


with  a  zero  eigen- 


The  vector  V'i(^)  [Hi(^)]  is  the  eigenvector  of  L{k,k3{k))  ^3(^)) 

value,  the  vector  V'2(^)  [^2(^)]  is  the  eigenvector  of  Z(fc, -A:3(^))  [X(^,  — fe3(^))]  with  a  zero 
eigenvalue  and  the  vectors  ^3, 4, 5(^)  [^3,4,5(^)]  are  the  eigenvectors  of  Z(0,  A:2,  fcs)  [X(0,A:2,^3) 

with  a  zero  eigenvalue. 


The  eigenvectors  which  correspond  to  kj  =  0,  (V3,4,5(^)  and  l=,3_4,5(X)),  represent  waves 
coming  from  —  oo  to  oo,  in  the  half  space  coordinates,  parallel  to  the  surface  dR^  in  the 
analysis  problem,  and  from  oo  to  -oo  in  the  adjoint  problem  and  therefore  do  not  play  a 
role  in  the  analysis. 

The  eigenvectors  which  correspond  to  k^  =  k3{k)  and  ks  =  k3(k)  are  relating  changes  in 

the  designed  surface  with  the  flow  field.  However,  only  Vi(k)  and  E:i(A:)  are  consistent  wdth 
the  far-field  boundary  conditions  as  discussed  in  Sec  4.3. 

We  look  for  amplitudes  of  the  vector  solutions  X’(^)  in  (5.10)  and  A{k)  in  (5.11)  consistent 
with  the  boundary  conditions.  From  the  boundary  condition  of  the  state  equation  (5.3)  we 
get  (see  (3.12)  together  with  Fi  in  (5.15)) 


kz^{k)  =  uikia{k). 


(5.17) 


The  boundary  condition  of  the  costate  equation  (5.8)  (see  (3.17)  together  with  the  definition 
of  5i  in  (5.16))  implies 

(?k3^{k)  —  —2kipu^{k)  (5.18) 

and  the  design  equation  (5.4)  implies 

g{k)  =  pv^iklc^'^ik)  +  l.o.t.  (5.19) 


Finally,  substituting  of  (5.17)  in  (5.18)  and  of  (5.18)  in  (5.19)  yields 

kt 


H{k)  =  -2/u 


2, .4  H 


k^k 


2p'^u^ 


31^3 


\{l  -  M^)kj  +  kl\' 


(5.20) 


6  Analysis  of  the  Hessian 

In  the  previous  sections  the  symbol  for  an  approximate  Hessian  near  the  minimum  was 
obtained,  namelv, 

H{kuk2)  =  2p^u%{ki,k2)  (6.1) 

j^4 

Hkuk2)  =  1^2(1 

The  fact  that  the  same  Hessian  is  obtained  both  for  the  full  potential  equation  and  for  Euler 
equations  implies  that  for  the  purpose  of  developing  new  optimization  algorithms  it  is  enough 
to  consider  the  full  potential  equations.  Although  the  Euler  equations  presents  additional 
difficulties  compared  with  the  full  potential  equation  for  the  analysis  problem  no  additional 
ones  exist  as  far  as  the  optimization  is  concerned. 
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The  Hessian’s  symbol  in  the  discrete  space  can  be  obtained  in  an  analog  way.  The  result 
will  then  depend  on  the  specific  discretization  we  use  to  solve  the  equations.  Gne  can  get 
a  crude  approximation  of  the  discrete  Hessian  by  replacing  the  wave  number  k\  with  the 
discrete  wave  number  di  (and  k2  with  O2), 


(h  h\ 

{ki,k2)  =  {-r;-r) 

h\  1I2 


where  {hi,h2)  are  the  mesh-sizes  in  the  {x,y)  directions  respectively. 


6.1  Two  Dimensions 

In  two  dimensions  the  X2  direction  does  not  exist  and  thus  in  the  Fourier  space  we  set  k2  =  0 
in  Eq.(6.2)  resulting  in 

Hh)  =  (6.3) 

Substitution  of  in  (6.3)  implies  that  the  condition  number  of  the  Hessian  scales  with 

the  grid  mesh-size  as  0{-^).  Thus,  the  Hessian  is  ill-conditioned  and  its  condition  number 
increases  quadratically  in  the  discretization  parameter  as  the  grid  is  refined.  Therefore,  a 
good  estimate  of  the  Hessian  is  required  to  obtain  a  fast  convergence  of  the  optimization 
process. 


6.2  Three  Dimensions 

In  three  dimensions  the  properties  of  the  Hessian  are  determined  by 

^  \kj{l  -  J^P)  +  kj\' 

For  fixed  wave  number  in  the  stream  direction,  ki  —  const,  the  symbol  approaches  zero 
as  the  wave  number  k2  approaches  infinity.  This  means  that  the  cost  function  is  nearly 
flat  with  respect  to  perturbations  in  the  shape  which  are  highly  oscillatory  in  the  direction 
perpendicular  to  the  flow.  This  might  explain  recent  numerical  results  showing  that  the 
wing  surface  is  likely  to  develop  oscillations  in  the  spanwise  direction  [9]  .  Note  that  the 
above  oscillations  do  not  appear  in  the  2D  problem  which  indicates  that  the  3D  aerodynamic 
optimal  shape  design  is  inherently  a  more  difficult  problem  than  the  2D  problem. 

7  Regularization 

As  discussed  in  Sec. 6. 2  the  three  dimensional  problem  is  ill-posed  and  as  a  result  oscillations 
are  expected  to  appear  in  the  direction  perpendicular  to  that  of  the  flow.  In  order  to 
eliminate  this  phenomenon  it  is  necessary  either  to  penalize  the  cost  function  or  alternatively 
to  represent  the  shape  as  a  finite  sum  of  base  functions  which  are  smooth  in  the  direction 
perpendicular  to  that  of  the  flow  (roughly  the  spanwise  direction). 
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In  this  section  we  propose  a  natural  penalty  of  the  cost  function  which  will  prevent  the 
oscillations  without  increasing  the  computational  cost  of  the  preconditioner.  The  penalty 
term  has  the  following  form, 

Fr  =  rj  l^iTafda  (7.1) 

where  77  is  a  positive  parameter  and  T  is  an  operator  whose  symbol  is 


fih,k2)  = 


kl 


\kl{l-M^)  + 


This  results  in  a  modified  Hessian  given  by 

+  (7.2, 

which  reaches  asymptotically  a  constant  as  k2  increases,  for  fixed  value  of  ki .  Thus,  for  large 
values  of  |fci|  +  |it2|  vve  have  \H{ki,k2)\  >  ^  >  0  for  some  positive  6,  which  implies  that 
for  the  modified  problem  the  high  frequency  perturbations  in  the  shape  are  well  behaved. 
As  the  shape  T  converges  to  the  minimum,  t]  can  be  decreased,  resulting  in  a  solution 
which  is  regularized  in  the  direction  perpendicular  to  the  flow.  The  implementation  of  the 
regularization  is  discussed  in  the  next  section. 


8  Preconditioning 

In  Sec.6  we  concluded  that  aerodynamic  optimization  are  ill-conditioned  problems.  Therefore 
having  second  order  information  is  crucial  for  effective  convergence.  However,  the  explicit 
solution  of  a  Newton  step,  Ha  =  —g,  requires  to  compute  the  Hessian,  H,  explicitly  and  then 
to  invert  it.  This  will  become  practically  impossible  for  a  realistic  aerodynamic  optimization 
problem  computed  numerically  on  a  fine  mesh  and  having  a  large  number  of  design  variables. 
Using  low  rank  quasi-Newton  methods,  such  as  BFGS,  will  deteriorate  as  the  number  of 
design  variables  increases.  We  therefore  suggest  to  approximate  the  Newton  step  in  the 
differential  level  (infinite  dimension),  using  the  Hessian’s  symbol,  and  then  to  project  the 
result  onto  the  finite  dimensional  design  space  which  is  used  in  practice. 

The  equation  defining  the  Newton  direction,  for  a,  in  Fourier  space  is  given  by 

H{ki,  k2)a{ki,  k2)  =  -g{ki,  k2) 

or  explicitly  by  using  the  symbol  (6.1) 

{2p\^kt  +  2gkl)a{k,M)  =  -|^i(l  -  M^)  +  kl\g{k^,k2)-  (8.1) 

The  symbol  in  the  right  hand  side  of  equation  (8.1)  corresponds  to  a  non-local  operator  in 
the  real  space.  The  term  which  multiplies  7/  accounts  for  the  regularization  penalty  term 
discussed  in  the  previous  section  (77  should  be  set  to  zero  in  two  dimensional  problems). 
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8.1  Preconditioner  for  the  Small  Disturbance  Problem 

Using  the  relation  between  ^3,  ki  and  A:2, 


we  arrive  at 

{2p^u^k^  +  2r)kl)a{ki,  k2)  =  {ikz){ikz)g{ki,  k2).  (8.2) 

This  implies  the  following  equation  in  real  space, 

(8.3) 

The  terms  which  multiply  fx  are  added  both  to  ensure  a  unique  solution  to  Eq.(8.3)  and 
to  account  for  the  low-frequencies.  Note  that  Eq.(8.2)  is  a  good  approximation  for  the 
symbol  of  the  Newton  equation  only  in  the  high-frequencies.  In  the  low-frequencies  the 
terms  multiplying  //,  in  Eq.(8.3),  are  dominant  and  result  in  a  steepest  descent  step,  while  in 
the  high-frequency  regime  they  diminish  and  a  Quasi-Newton  step  is  taken.  The  term 
satisfies  the  following  coupled  PDE  system 


on  1R\ 

(8.4) 

OZ 

on  d]R\ 

(8..5) 

on  1R\ 

(8.6) 

li 

on  dJR\. 

(8.7) 

We  also  require  that  the  solution  be  bounded  and  that  in  the  supersonic  regime  {M  >  1)  no 
waves  propagate  in  the  positive  x-direction  in  Eq.(8.5)  and  that  no  waves  propagate  in  the 
negative  x-direction  in  Eq.(8.7).  Note  that  the  operator  T  in  Eq.(7.1)  need  not  be  explicitly 
evaluated.  By  adding  the  rj  term  in  (8.2)  and  solving  Eqs.(8.3-8.7)  we  account  for  such  an 
operator.  A  similar  preconditioner  can  be  derived  for  the  small  disturbance  Euler  equations. 


8.1.1  Purely  Subsonic  Flow 

In  a  purely  subsonic  flow  it  is  unnecessary  to  go  through  the  above  procedure  since 

\kl{l  -  M^)  +  kl\  =  kl{l  -  M^)  +  kl 


Therefore  the  following  PDE  should  be  solved  on  T  only: 

,<9^(5;  ^  d'^a 


+  ga  =  {l-  M^)^^  +  ^  + 


dx^ 


d^y 


d^g  d^g 

d^x  d^y 


(8.8) 
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8.2  Practical  Implementation 

In  practice  equations  (8.4-8. 7)  are  replaced  by  full  potential  equations  defined  on  the  domain 
O  as  shown  in  the  following.  Let  be  the  solution  of  the  full  potential  equation  (4.2).  We 
define 

^(^1  =  tp)  -\- 

w'here  £  <  1.  We  claim  that  the  preconditioning  equation  for  d  is  given  by 

4  2  1  (9 

-  2?/(e  •  -f -  V?) -f  on  F,  (8.9) 

where  e  is  a  unit  vector  perpendicular  to  V</?  and  to  n  (roughly  in  the  spanwise  direction); 
X  e  =  n.  The  function  is  the  solution  of  the  following  coupled  PDE  system: 

=  0  on  0 

on  r  (8.10) 

on  dfl  —  r 

Vp(^(2))y^(2)  ^0  on  D 

-(- (^  on  r  (8-11) 

^(2)  _  dfloj 

where  dO-of  stands  for  the  far  field  outflow  boundary.  The  solution,  d,  of  Eq.(8. 9-8.11)  is  the 
preconditioned  search  direction  to  be  used  in  optimization  algorithms.  This  search  direction 
will  avoid  oscillations  in  the  shape  in  the  direction  perpendicular  to  the  flow  and  will  require 
many  fewer  optimization  steps  to  solve  the  problem.  For  the  Euler  equations  an  analog  of 
Eqs.(8.10-8.11)  can  be  derived  using  the  corresponding  small  disturbance  preconditioner. 

8.3  Implementation  in  a  Finite  Dimensional  Design  Space 

In  a  finite  dimensional  subspace  we  replace  d  in  Eq.(8.9)  by  (see  (2.15)) 

d(f)  =  ^ajfj{x). 
i=i 

We  then  take  the  inner  product  of  Eq.(8.9)  with  fk  iov  k  =  1  ■  •  ■  N ,  resulting  in  a  linear  set 
of  equations  for  aj  : 

{h  +  nl)a  =  -6  -f  ng,  (8.12) 

where  I  is  the  unit  matrix  and  yu  is  a  positive  parameter,  as  in  Eq.(8.3).  h  is  &n  N  x  N 
matrix  in  which  an  element  hij  is  given  by 

=  (2p2(v^  .  v)Vt  -  2r/(e  •  V)  , 
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Q  is  the  finite  dimensional  design  variable  a  =  (ai,  •  •  • ,  Qjv),  and  an  element  in  the  right 
hand  side  vector,  6,  is  given  by 


k  =  — 
e 


dn 


(1^(1)  _  ^  f. 


L2(r) 


The  solution  of  Eq.(8.12)  is  the  preconditioned  search  direction  replacing  •  •  •  ,5'Ar). 


9  Conclusions 

In  this  work  new  local  mode  analysis  for  optimal  shape  design  problems  which  are  governed 
by  PDEs  was  developed.  The  analysis  was  applied  to  aerodynamic  shape  optimization 
problems  governed  by  the  full  potential  and  the  Euler  equations.  The  analysis  was  done  in  the 
infinite  dimensional  design  space  where  arbitrary  changes  in  the  wing’s  shape,  in  the  normal 
direction,  are  allowed  during  the  course  of  optimization.  In  this  case  the  Hessian  is  an  infinite 
dimensional  operator  defined  on  a  space  of  functions  and  its  eigenvalue  distribution  served  to 
study  the  well-posedness  of  the  optimization  problem  as  well  as  for  deriving  preconditioners 
to  accelerate  the  numerical  convergence  of  gradient  based  methods.  In  practice,  however,  a 
finite  dimensional  design  space  is  commonly  used  for  which  the  Hessian  is  a  finite  dimensional 
matrix.  The  application  of  the  infinite  dimensional  analysis  to  finite  dimension  was  obtained 
by  a  simple  projection. 

The  analysis  is  local  and  uses  freezing  the  coefficient  to  obtain  a  problem  in  half  space  with 
constant  coefficients,  where  Fourier  techniques  are  employed.  The  eigenvalue  distribution 
of  the  Hessian  is  analyzed  by  computing  its  Fourier  symbol.  It  was  shown  that  for  two 
dimensional  flow  the  Hessian  is  a  second  order  differential  operator  defined  on  the  designed 
boundary.  In  three  dimension  the  Hessian  is  a  pseudodifferential  operator  (non-local)  and  its 
properties  are  much  more  complex.  For  both  the  full  potential  and  Euler  flow  the  symbols 
of  the  Hessian  are  identical.  Therefore,  the  complexity  of  the  optimization  problems  is  the 
same  for  both,  although  the  analysis  problem  for  the  Euler  equations  is  more  difficult. 

The  symbol  of  the  Hessian  implies  that  the  three  dimensional  problems  are  ill-posed,  and 
arbitrary  oscillations  in  the  shape  can  develop  in  the  direction  perpendicular  to  that  of  the 
flow  (roughly  the  spanwise  direction).  This  explains  recent  numerical  results  showing  that 
the  wing  surface  is  highly  oscillatory  in  the  spanwise  direction  [9]  .  A  regularization  that 
involves  smoothing  only  in  that  direction  was  introduced  and  analyzed.  Also  note  that  if 
the  problem  were  to  minimize  drag,  rather  then  matching  the  pressure  distribution,  then 
oscillations  in  the  spanwise  direction  are  not  likely  to  appear  since  they  will  increase  the 
surface  area  of  the  wing  resulting  in  an  increase  in  the  drag. 

The  explicit  form  of  the  symbol  of  the  Hessian  also  implies  that  these  minimization 
problems  are  ill-conditioned  and  their  condition  number  increases  quadratically  with  the 
dimension  of  the  design  space.  Therefore  gradient  descent  method  will  be  inefficient  and 
second-order  information,  by  approximating  the  Hessian  (or  its  inverse),  is  essential  for  fast 
convergence.  However,  low  rank  quasi-Newton  methods,  such  as  BFGS,  will  deteriorate 
as  the  number  of  design  variables  increases.  New  preconditioners  which  approximate  the 
inverse  of  the  Hessian  are  proposed.  Their  numerical  implementation  will  be  presented 
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elsewhere.  These  preconditioners  are  of  low  computational  complexity  for  two-dimensional 
flow  and  for  purely  subsonic  flows  in  three  dimensions.  In  non-subsonic  three-dimensional 
flow  the  preconditioning  involves  the  solution  of  the  full  potential  equation  twice  per  each 
optimization  step,  though  a  substantial  decrease  in  the  number  of  optimization  steps  required 
to  reach  the  minimum  is  anticipated.  The  preconditioning  of  the  suggested  regularization 
is  straightforward  and  requires  negligible  additional  computational  work.  In  case  a  flnite 
dimensional  design  space  is  used  the  preconditioning  requires  the  solution  of  a  linear  set  of 
equations. 

Applications  of  similar  analysis  to  aeroelastic  optimization  are  discussed  in  [18]. 

Acknowledgment 

The  authors  would  like  to  thank  R.  M.  Lewis  for  his  important  comments  on  the  manuscript. 


References 

[1]  0.  Pironneau.  Optimal  Shape  Design  for  Elliptic  Systems,  Springer  Series  in  Compu¬ 
tational  Physics  (1983). 

[2]  B.  Mantel,  J.  Periaux  and  B.  Stoufflet.  AGARD-FDP-VKI  Special  Course,  Optimum 
Design  Methods  in  Aerodynamics,  April  25-29  (1994). 

[3]  A.  Jameson.  Aerodynamic  Design  Via  Control  Theory,  Journal  of  Scientific  Computing. 
3:233-260  (1988). 

[4]  A.  lollo  and  M.  Salas.  Contribution  to  the  Optimal  Shape  Design  of  Two-Dimensional 
Internal  Flows  with  Embedded  Shocks,  ICASE  Report  No.  95-20  (1995). 

[5]  S.  Ta’asan,  G.  Kuruvila  and  M.  D.  Salas.  Aerodynamic  Design  and  Optimization  in 
One  Shot,  30th  Aerospace  Sciences  Meeting  &:  Exhibit,  AIAA  92-0025,  Jan.  (1992). 

[6]  S.  Ta’asan.  Trends  in  Aerodynamics  Design  and  optimization:  A  Mathematical  View¬ 
point.  AIAA  95-1730,  12th  AIAA  Computational  Fluid  Dynamics  Conference,  July 
19-22  (1995). 

[7]  A.  Dervieux,  J.  Male,  N.  Macro,  J.  Periaux,  B.  Stoufflet  and  H.Q.  Chen.  Some  Recent 
Advances  in  Optimal  Shape  Design  For  Aeronautical  Flows,  Proceedings  of  “ECCO- 
MAS,  2nd  Computational  Fluid  Dynamics  Conference”,  Sep.  5-8  (1994). 

[8]  P.  D.  Frank  and  G.  R.  Shubin.  A  comparison  of  optimization-based  approaches  for  a 
model  computational  aerodynamics  design  problem.  J.  Comput.  Phys.  98,  74  (1992). 

[9]  W.  H.  Jou,  W.  P.  Huffman,  D.  P.  Young,  R.  G.  Melvin,  M.  B.  Bieterman,  C.  L.  Hilmes 
and  F.  T.  Johnson  (The  Boeing  Company).  Practical  Considerations  in  Aerodynamic 
Design  Optimization,  AIAA  95-1730,  12th  AIAA  Computational  Fluid  Dynamics  Con¬ 
ference,  July  19-22  (1995). 

[10]  J.  Hadamard.  Lessons  on  the  Calculus  of  Variation  (in  French),  Gauthier- Villards, 
Paris  (1910). 


19 


[11]  M.  J.  Lighthill.  A  new  method  of  two  dimensional  aerodynamic  design.  R  &  M  1111, 
Aeronautical  Research  Council  (1945). 

[12]  P.  E.  Gill,  W.  Murray  and  M.  H.  Wright.  Practical  Optimization,  Academic  Press 
(1981). 

[13]  E.  Arian.  Multigrid  Methods  for  Optimal  Shape  Design  Governed  by  Elliptic  Systems, 
Ph.D.  Thesis,  The  Weizmann  Institute  of  Science,  Israel  (1994). 

[14]  E.  Arian  and  S.  Ta’asan.  Multigrid  One  Shot  Methods  for  Optimal  Design  Problems: 
Infinite-Dimensional  Control,  ICASE  Report  No.  94-52  (1994). 

[15]  E.  Arian  and  S.  Ta’asan.  Shape  Optimization  in  One  Shot,  Optimal  Design  and  Con¬ 
trol,  Edited  by  .J.  Boggaard,  .J.  Burkardt,  M.  Gunzburger  and  .J.  Peterson,  Birkhauser 
Boston  Inc.  (1995). 

[16]  E.  Arian  and  S.  Ta’asan.  Smoothers  for  Optimization  Problems,  Proceedings  of  the 
Seventh  Copper  Mountain  Conference  on  Multigrid  Methods,  April  2-7  (1995). 

[17]  C.  Hirsch.  Numerical  Computation  of  Internal  and  External  Flows,  Vol.  2,  John  Wiley 
k  Sons  (1988). 

[18]  E.  Arian.  Analysis  of  the  Hessian  for  Aeroelastic  Optimization,  ICASE  Report  No. 
95-84  (1995). 


20 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


D  .  _ for  thi^  rolleciion  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 

Public  rep  £,;n4,;ninffthe  data  needed  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 


1.  AGENCY  USE  OUVf (Leave  blank)  2.  REPORT  DATE 

April  1996 _ 


4.  TITLE  AND  SUBTITLE 

ANALYSIS  OF  THE  HESSIAN  FOR  AERODYNAMIC 
OPTIMIZATION;  INVISCID  FLOW 


3.  REPORT  TYPE  AND  DATES  COVERED 

Contractor  Report  _  ' 


5.  FUNDING  NUMBERS 

C  NASl-19480 
WU  505-90-52-01 


6.  AUTHOR(S) 

Eyal  Arian 
Shlomo  Ta’asan 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science  and  Engineering 
Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  96-28 


9,  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23681-0001 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

NASA  CR-198328 
ICASE  Report  No.  96-28 


11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  BushneU 
Final  Report 

To  be  submitted  to  the  Journal  of  Computational  Physics. 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

U  nclassified-Unlimited 


1  12b.  DISTRIBUTION  CODE 


Subject  Category  64 

13.  ABSTRACT  (Maximum  200  words)  ...  ,  ,  j  r  xi.  r  ii  +  l  ^ 

In  this  paper  we  analyze  inviscid  aerodynamic  shape  optimization  problems  governed  by  the  full  potential  and 

the  Euler  equations  in  two  and  three  dimensions.  The  analysis  indicates  that  minimization  of  pressure  dependent 
cost  functions  results  in  Hessians  whose  eigenvalue  distributions  are  identical  for  the  fuU  potential  and  the  Euler 
equations.  However,  the  optimization  problems  in  two  and  three  dimensions  are  inherently  different.  While  the  two 
dimensional  optimization  problems  are  well-posed,  the  three  dimensional  ones  are  ill-posed.  Oscillations  in  the  shape 
up  to  the  smallest  scale  allowed  by  the  design  space  can  develop  in  the  direction  perpendicular  to  the  flow,  implying 
that  a  regularization  is  required.  A  natural  choice  of  such  a  regularization  is  derived.  The  analysis  also  gives  an 
estimate  of  the  Hessian’s  condition  number  which  implies  that  the  problems  at  hand  are  ill-conditioned.  Infinite 
dimensional  approximations  for  the  Hessians  are  constructed  and  preconditioners  for  gradient  based  methods  are 
derived  from  these  approximate  Hessians. 


NSN  7540-01-280-5500 


15.  NUMBER  OF  PAGES 

22 


16.  PRICE  CODE 

A03  _ 


18  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION 
OF  THIS  PAGE  OF  ABSTRACT  OF  ABSTRACT 

Unclassified  _ 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


14.  SUBJECT  TERMS  t.  i  t. 

Aerodynamic  Optimization;  Optimal  Shape;  Full  Potential  Equations;  Euler  Equations 

Hessian;  Preconditioning;  Regularization 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


